Self-Sustaining Oscillations in Complex Networks of Excitable Elements 
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Abstract 

Random networks of symmetrically coupled, excitable elements can self-organize into coherently oscillating states if the 
networks contain loops (indeed loops are abundant in random networks) and if the initial conditions are sufficiently random. 
In the oscillating state, signals propagate in a single direction and one or a few network loops are selected as driving loops in 
which the excitation circulates periodically. We analyze the mechanism, describe the oscillating states, identify the pacemaker 
loops and explain key features of their distribution. This mechanism may play a role in epileptic seizures. 

PACS numbers: 89. 75. He, 05.65.+b, 87.19.1j, 87.19.xm 



The coherent oscillation (CO) of a collection of units 
that are non-oscillatory on their own is relevant to 
biological and physical sciences: CO has been identi- 
fied and analyzed in populations of excitable biologi- 
cal cells (yeast j3 pancreatic cellsH, Dictyostelium 
discoideum[l[ and cultured heart cellsQ) and of excitable 
catalytic particles [5] l 6]|7]. In contrast with well-studied 
synchronization phenomena of self-oscillating units Q, in 
these cases the ability to oscillate derives from the in- 
teractions of the elements. CO can occur also on com- 
plex networks if the nodes are excitable or even 
monostableQ, provided that the network contains loops, 
and that the directional symmetry of couplings is some- 
how broken to allow a signal to propagate in a sin- 
gle direction around a loopjioj. Networks of these 
types include some neural [l2j and genetic regulatory 
networks. @ Some studies of excitable networks have 
been inspired by target and spiral waves in continuous 
media. 

In com plex networks, loops are both generic and 
abundant. |15j|[ll| While short loops of length L <C N 
(where N is the network size) are rare in large random 
networks, ones with L > In TV occur generically in num- 
bers growing exponentially with N. Their number also 
grows roughly exponentially with L up to a maximum at 
a length L m ~ N. (UPl 

In this letter we show that random excitable networks 
readily self-organize into a CO state following a transient 
phase during which one or a few loops are dynamically 
selected as driving (or pacemaker) loops. We describe 
the mechanism of signal propagation and gain an under- 
standing of the resulting distributions of the oscillating 
states and associated driving loops. 

We consider networks of diffusively coupled, excitable 
elements with dynamics described by the Bar model [l6| 
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N is the number of nodes, Ui and Vi are dynamical vari- 
ables, D is the coupling strength, and A t j is the adja- 
cency matrix. a, 6, and e are parameters, for which 
we adopt the values a — 0.84, b — 0.07, and e = 0.04. 
In the absence of coupling, the individual nodes dis- 
play excitable dynamics, with a stable equilibrium at 
(u, v) = (0,0) and an excitation threshold u t h ~ 0.1. 
The dynamical equations were integrated numerically us- 
ing a fourth-order Runge-Kutta algorithm with time step 
At = 0.1. 

For the topology, we chose the undirected random reg- 
ular network of degree 3 (RRN3), where nodes are ran- 
domly and symmetrically connected with the constraint 
that all have the same degree k = 3. The network size 
was N — 200. Other topologies and other values and 
distributions of k will be considered elsewhere. 

To excite the network, we used random initial condi- 
tions in which each node was either displaced from equi- 
librium with probability p, or left at (0, 0) with probabil- 
ity 1— p. For the nodes that were displaced, initial values 
of u and v were distributed randomly and independently 
within the intervals 0.2 < u < 0.9, < v < 1, thus plac- 
ing them above the excitation threshold but with some 
phase randomness. (This randomness proves important, 
as discussed below.) After determining that the results 
were largely insensitve to the value of p, we subsequently 
took p — 0.5. 

Integrating the equations of motion (JTJ) , we found two 
possible outcomes: either the system relaxed rapidly to 
the quiescent state with all nodes at the fixed point or it 
reached a coherently oscillating state (COS) like the one 
illustrated in figHJ In the COS, all nodes fire at the same 
frequency, with a fixed phase relationship among them. 
To examine the phase relations, we define a firing time U 
of the zth node as the time (interpolated linearly between 
discrete time steps) at which Ui crosses from below 0.5 to 
above 0.5. The interval between successive firing times 
(interspike interval or ISI) converges to a stable common 
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FIG. 1: (A) Probability that a random initial condition with 
excitation probability p = 0.5 relaxes to an oscillatory state, 
averaged over 1000 initial conditions of 5 network realizations. 
(B) Scatter plot of oscillation periods T. Periods are longer 
and more variable at weaker coupling. 



value for all nodes, which is the oscillation period T. We 
considered the system to have converged to a COS if the 
standard deviation of the ISI's over the network remained 
below 10 -4 for more than 100 time units. As shown in 
fig HI COS arc highly probable for a range of coupling 
strengths 0.1 ;$ D ;$ 0.6, while the periods decrease with 
increasing D. 

The local structure of the oscillatory state is illustrated 
by fig[2}3, which shows the firing pattern of an arbitrarily 
chosen node and its three neighbors. As each node 
fires once per oscillation, one can define a time delay Ay 
(with -T/2 < Aij < T/2) for each link in the network 
as the difference of firing times t% — tj between the nodes 
at its two ends, modulo the oscillation period T. With 
the establishement of stable signed delays, the initially 
undirected network self-organizes into a directed one in 
which the signal propagates in only one direction along 
each link. Since each node i must be excited by one 
of its neighbors, it must have at least one incoming link 
(positive Ay ) but can have between and 2 outgoing 
links (negative Ay). The node illustrated in figEP is 
a "diverging" node with one incoming and two outgoing 
links. Since the total numbers of incoming and outgoing 
links in the network must balance, all nodes cannot be 
converging. This implies that CO can only occur if a 
single firing neighbor suffices to excite a node. In this 
case, the firing of a node at the "upstream" end of a link 
guarantees that the one at the " downstream" end will fire 
within a certain time period, provided the downstream 
node is not refractory when it receives the input. If the 
downstream node receives a second input before it fires, 
it will be pulled over the threshold more quickly and fire 
sooner. Converging inputs thus account for most of the 
variation in transmission delays. 

Liao et al. [H suggested a way of identifying the main 
transmission pathways by selecting the so-called domi- 



(B) 



k 'J 


; I ; I 






If; p 




ft ; 






v •. 






'Ml! : 






:/ 


.1 


7 


,/ 


Wi 


l\/ 


- 




FIG. 2: (color online): Example of a periodic network oscil- 
lation for coupling strength D — 0.2. (A) Raster plot of all 
nodes' activity. Gray level is proportional to Ui(t). The 
network quickly settles to a steady oscillation in which all 
nodes fire with the same period in a fixed phase relationship 
to each other. (B) u(t) for an arbitrarily chosen node (solid 
line) and its three neighbors (incoming: dotted and outgoing: 
dashed lines). 
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FIG. 3: Branched circle structure of the pruned network. 
Following dominant phase-advanced driving links backwards 
from any node inevitably leads to a unidirectional loop from 
which the signal radiates along branches. 



nant phase-advanced driving (DPAD) links, which then 
leads to identifying the pacemaker or driving loops (DL). 
The DPAD for each node is the incoming link with the 
largest delay. A justification for considering the earli- 
est input to be the most important is that one input is 
sufficient to guarantee the node's firing, and additional 
inputs can only affect the timing Pruning the network 
to include only DPAD links simplifies it to a so-called 
branched circle structure |9| consisting of trees attached 
to unidirectional DLs as illustrated in fig. [3J 

To study the statistics of the COS's and their basins 
of attraction we integrated eqs. (HJ with D = 0.11 for 
1000 different initial conditions (200 each for five different 
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RRN3 realizations). All but one of these initial condi- 
tions converged to a COS. We then measured the fir- 
ing delays and pruned the network as described above to 
identify the driving loops. We found up to four distinct 
DLs in each COS, all mutually entrained to oscillate with 
the same period. 50.9 percent of the COS had only one 
DL, with the probability of more loops decreasing mono- 
tonically with number. In each case, we measured the 
length L of the shortest DL. The distribution of these 
lengths is shown in figure HJ\. As figure [4j3 shows, the 
oscillation period T is correlated with L, but the data fall 
into clusters separated according to the number of pulses 
circulating simultaneously around the loop, which we call 
the multiplicity M. M can be measured by adding up 
the transmission delays along the driving loop to get the 
time for a single pulse to make one circuit, and dividing 
this by the oscillation period. FigSfU shows that the 
data collapse into a single band when one plots T vs. 
L/M. One can interpret the slope of this band as the 
typical delay added by a link. A notable feature of this 
plot is that the best fit line does not pass through the ori- 
gin as one would expect if waves of excitation travelled at 
a constant velocity independent of the loop length. This 
feature becomes even more noticeable when we examine 
data at stronger coupling (where transmission delays are 
shorter). Data for D = 0.6 look qualitatively similar to 
those shown in figure 21 except that the minimum value 
of L is 6 rather than 4, and the slope of the plot of T 
vs. L/M is considerably smaller, as expected, while the 
intercept is only slightly smaller. The resulting overall 
spread of oscillation periods is correspondingly less. At 
larger coupling strength, the transmission delay is smaller 
and evidently plays less of a role in setting the period. 

Several features of our results merit discussion. First, 
we note the ease with which the network is induced to 
oscillate: the probability of obtaining oscillations from 
random initial conditions is nearly 1 within a range of 
coupling strengths, and there is a large number of distinct 
oscillatory attractors as attested by the scatter of periods 
and the variety of DLs. Second, no external driving is 
needed to sustain CO activity, unlike the case considered 
by [IlJ- Unlike cases where oscillation resulted from 
adding directed shortcuts to a spatial network (7j [l3T| fl7| . 
here the connections are symmetric. This symmetry 
is only broken dynamically, producing different directed 
networks and selecting different driving loops depending 
on the initial conditions. 

A non-trivial distribution of the initial conditions in 
phase space is evidently crucial: initial conditions where 
one or more nodes were excited synchronously (displaced 
to the same point above the excitation threshold) failed 
to produce oscillations. A single excited node on a loop 
produces two wavefronts of excitation propagating in op- 
posite directions, which will annihilate elsewhere on the 
loop. If an excited node adjoins a refractory one, how- 
ever, the signal is blocked from moving in one direction 
and a unidirectional loop can be established. [10] Evi- 
dently the random initial distribution we have used is 
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FIG. 4: (color online): Pacemaker loops and oscillation pe- 
riods for D = 0.11, for 1000 different initial conditions of 5 
network realizations. (Features are qualitatively similar for 
COS at D — 0.6) (A) Histogram of lengths of shortest driv- 
ing loops. (The shape is qualitatively similar if all loops are 
included.) A semilogarithmic plot of the cumulative distribu- 
tion (not shown here) reveals an exponentially decaying trend. 
(B) Oscillation period T vs. shortest loop length. Different 
symbols show loops of different multiplicities M. Note that 
the clusters at different M values correspond to the local peaks 
in panel (A). (C) The data from (B) collapse into a single 
band when one plots T vs. L/M. The best fit line shown here 
is T = 1.21(L/M) + 5.26. For D = 0.6, the corresponding 
best fit is T = 0.25(L/M) + 2.85. 



sufficient to allow this. 

As shown in figure [1] the probability of obtaining os- 
cillations from a random initial condition jumps rapidly 
from to 1 at a lower threshold coupling strength Di = 
0.1, and falls off more gradually for D ^ 0.6. The lower 
threshold represents the minimum coupling necessary for 
a node to be excited by the firing of a single neighbor. 
The slower decay of the oscillation probability at D Si 0.6 
can be explained by the shortening of transmission delay 
compared to the refractory period. In the cases where os- 
cillation fails, it does so due to an avalanche that spreads 
so rapidly through the network that almost all nodes are 
firing at once. The activity burns out quickly when there 
are no nodes remaining that are not already firing or re- 
fractory. 

The oscillation period is bounded from below by the 
refractory period. This also implies a lower bound on 
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the length of a DL depending on the maximum delay per 
link: oscillation cannot be sustained if a pulse travelling 
along a loop returns to a given node while that node is 
still refractory. The observed lower bound on the loop 
size, -L m in = 4 in the case D = 0.11, in fact increases to 
6 when D is increased to 0.6, because increasing D de- 
creases the transmission delay. More generally, the above 
argument explains the observed lower bound on the ratio 
L/M, which represents the spacing between pulses on a 
DL. L/M, on the other hand, is also apparently bounded 
from above, as evident from figure |3p. A possible expla- 
nation for this lies in the transient process by which a 
particular loop becomes established as the primary DL. 
During this transient phase, many pulses are propagat- 
ing in an unsynchronizcd manner through the network, 
and a number of potential pacemakers are competing. A 
large gap between circulating pulses is therefore likely to 
be filled in. This is analogous to oscillations in excitable 
spatial media, where pacemakers with the shortest peri- 
ods dominate. 

The distribution of driving loops (figure decays 
with L, despite the fact that the number of loops present 
in the network grows exponentially with L. There is ev- 
idently a strong bias in favor of shorter loops. This bias 
appears to be a statistical effect inherent in the topol- 
ogy of the directed branched circle network rather than 
strictly a dynamical effect. To check this hypothesis, we 
compared our results with a "null model" that generates 
branched circle networks independently of the oscillatory 
dynamics. For the null model, we generate directed net- 
works as follows. First, assign a random direction to 
each link in an RRN3, and then, where necessary, re- 
verse directions of some links to ensure that every node 
has at least one incoming link. For each node, we then 
randomly choose one among its incoming links to be the 
dominant one, and prune the others if the node is con- 
verging. The result is a network satisfying the same 
topological constraints as the pruned network of DPAD 
connections, but generated by a random process rather 
than a dynamical one. The distribution of the shortest 



loop length L for an ensemble of such networks decays 
exponentially as it does in the DPAD network. Simply 
by assigning directions and then pruning, one samples 
the loops of the original undirected network unevenly. 
Longer loops have more chances either to fail to be uni- 
directional when directions are assigned, or to be elimi- 
nated by pruning. 

The driving loops uncovered by the DPAD reduction 
explain part but not all of the variation in oscillation pe- 
riods (see figHJ). Secondary (non-dominant) connections 
alter the firing times, and are also responsible for the en- 
trainment of multiple DL when the latter occur. If the 
connections that we prune for the sake of analysis were 
truly removed from the network, it would be impossible 
for the resulting disconnected components to synchro- 
nize. It is also worth emphasizing that when M > 1, 
the multiple pulses are precisely evenly spaced as they 
travel around the loop, so that the interval between fir- 
ings is constant for a given node. This is another form 
of entrainmcnt that we would not expect if the pruned 
connections were truly absent. 

The dynamical organization of an undirected network 
giving rise to a directed one with considerably differ- 
ent loop statistics provides an intriguing case study of 
the distinction between "functional" networks defined 
by dynamical interactions and the underlying structural 
networks of hard- wired connections. Functional net- 
works of observed correlations of activity among brain 
regions have been studied as a tool to infer the underly- 
ing architecture [lSjj , but the phenomenon discussed here 
illustrates that functional and structural networks may 
have markedly different properties. 

Unlike the cases of interacting excitable cells where the 
oscillation is treated as a collective phenomenon in 
the case of a fixed network specific signalling pathways 
can be more readily identified. The networks studied 
here share some essential features with brains, and it is 
plausible that the mechanism described here plays a role 
in epileptic seizures. [l9j 
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